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Monte Carlo simulations are performed for hard ellipsoids for a number of values of its semi-axes 
in the range c/a G {0.1, 10}. The isotropic phase results are compared to the Vega equation of state 
[Mol. Phys. 92 651 (1997)]. The position of the isotropic-nematic transition is also evaluated. The 
biaxial phase is seen to form only after the previous formation of a discotic phase. 



INTRODUCTION 

One of the first ever models used in computer sim- 
ulation studies of the fluid phase was the hard sphere. 
The next most obvious choice of model is the hard ellip- 
soid; an affinc transformation of the hard sphere. The 
orientation of the model now becomes a variable within 
the simulation. Having a model that incorporates ori- 
entation facilitates the study of orientationally ordered 
phases, such as those associated with liquid crystals. 

As with hard spheres [l|, 0, 0, 0] the first simulations of 
hard ellipsoids were performed for two-dimensional sys- 
tems [1, 0] . Hard spheres are only capable of forming two 
phases; the fluid and the solid. There is no 'gas-liquid' 
transition, due to the lack of attractive forces. How- 
ever, the hard ellipsoid system also has a plastic crys- 
tal (for small axis ratios) and a nematic phase. Prolate 
and prolate-like biaxial models, i.e. models of the form 
a < b < c with b < ^fac where a, b and c arc the semi- 
axes of the ellipsoid, are able to form a uniaxial nematic 
phase, often denominated as N+. Oblate and oblate- 
like biaxial models, with a < b < c and b > \fac, can 
form a uniaxial 'discotic' phase (iV_). A tentative phase 
diagram for uniaxial ellipsoids was proposed by Frenkel 
and co-workers 0, S, 0, fiol. As well as these nematic 
and discotic phases, Preiser [11| predicted that 'long-flat' 
molecules could form a biaxial phase TVb, which has re- 
cently been discovered experimentally [l^, [l3| ■ 

The study of hard bodies provide reference systems for 
use in perturbation theories which add long range attrac- 
tive interactions IJ]. They can also form the monomer 
units of larger molecules [l^,[lB| • Hard ellipsoids continue 
to be of interest [13 1 having recently been shown that 
a maximally random jammed (MRJ) packing fraction of 
= 0.7707 is possible for models whose maximal aspect 
ratio is greater than \fZ. [11, [l^ , Such high packing frac- 
tions were obtained using an 'event-driven' molecular dy- 
namics code [2ll |22| and were confirmed experimentally 
using latex particles and the like 13, 2^ . 

Biaxial ellipsoids have, however, received compara- 
tively little attention, simulations having been performed 
principally by Allen 25 1 and by Camp and Allen [2^. Ex- 
tensive simulation results have been presented previously 
for uniaxial hard ellipsoids, notably those of Prenkel and 
Mulder for c/a in the range [1/3, 3] [^ In this publication 



a number of simulations are performed for biaxial ellip- 
soids within the region c/a for [0.1, 10]. Also in this work 
the isotropic equation of state is compared to theory. 



SIMULATION TECHNIQUE 

Standard Metropolis Monte Carlo sampling was used 
[l[. The decision to reject or accept a trial move for hard 
bodies is based simply on whether two bodies overlap 
or not. Por hard spheres (a = 6 = c) this criteria is 
trivial; if the distance between two bodies is less than 
twice the radius then they overlap. Por hard ellipsoids 
the situation is more complicated, having to take into 
account the orientations of the ellipsoids as well as the 
distance between them. One of the first overlap criteria 
was proposed for two-dimensional ellipsoids by Vieillard- 
Baron [1,[23]. Perram and Wertheim produced an overlap 
algorithm for three-dimensional ellipsoids [2^. In this 
study the Perram- Wertheim criteria was used. 

All of the Monte Carlo simulations were performed in 
the NpT ensemble (with ksT =1), having cubic bound- 
ary conditions, with the system comprising of N = 343 
hard ellipsoids. The runs consisted of between 75 and 
150 kilocycles for equilibration followed by a further 75- 
150 kilocycles for the production of thermodynamic data. 
Each simulation was initiated from the final configuration 
of the previous, lower pressure, run. Throughout the sim- 
ulations the uniaxial order parameter 5*2 was monitored 
for each of the three axis by calculating 



52-(i(3cos2 0,-i; 



(1) 



where 9i is the angle between ellipsoid i and the director 
vector [2i,[33,[3lj. 

The prolate parameters are typified by a = 6, the 
oblate parameters by 6 = c and the biaxial parameters 
by a 7^ 6 7^ c where in this work a < 5 < c (in this work 
we define a = 1). 

Boublik and Nezbeda have shown [s^] that hard con- 
vex bodies can be described using three geometric expres- 
sions; the volume, V, the surface area, S and the mean ra- 
dius of curvature, R. The second virial coefficient, B2 for 
any hard convex body is given by the Isihara-Hadwiger 
formula 33, 33|, and provides a measure of the excluded 
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volume of a 'molecule' due to the presence of a second 
molecule, 



B2=RS + V, 



B2 
V 



= 1 + 3a, 



if one defines the 'non- sphericity' parameter 35 1 

RS 



(2) 
(3) 

(4) 



Values for i?, S* and V for ellipsoids are provided in Ap- 
pendix A. 



RESULTS 
Isotropic phase equation of state 

Accurate equations of state (EOS) are necessary if the 
hard ellipsoid fluid is to be used as a reference system 
in perturbation theories. A number of equations of state 
have been proposed over the years to reproduce the be- 
haviour of the isotropic pha se of the ellipsoid system; no- 



tably those of Nezbeda [36| , Parsons [37[ , and Song and 
Mason jsl]. It should be noted that all of these equa- 
tions of state implicitly assume that the EOS is symmet- 
ric with respect to prolate/oblate models, since they are 
built either directly using B2, or indirectly using a. This 
is because i?2(l x 1 x x) = -62(1 x x x x). 

In order to reproduce simulation results, higher virial 
coefficients are required, where it has been shown that 
the prolate/oblate symmetry breaks down starting from 
-B3 [39I . kdl , as used in the uniaxial EOS of Maeso and 
Solana [4l[ or the biaxial EOS of Vega [i^l- In this work 
the simulation results are compared to the biaxial Vega 
equation of state. The Vega EOS is given by 



Z= 1 



4 I 

1 -4y- 



- 18.3648y^ 



28.2245?;^) , (5) 



where Z is the compressibility factor and y is the volume 
fraction, given hy y = pV where p is the number density. 
The virial coefficients are given by the fits 

El = 10 + 13.094756a' - 2. 073909r' + 4.096689a'2 
+ 2.325342t'2 - 5.791266aV, (6) 
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FIG. 1: Plot of the simulation results for the biaxial molecules 
1x2x4, 1x2x5, 1x2x6 and the isotropic state points 
for 1x2x8. The lines are the Vega equation of state. 



and 



b: = 



28.2245 + 21.288105a' + 4.525788r' + 36.032793a'2 

59.0098t'2 - 118.407497a'T' 

24. 164622a' V + 139.766174a'r'2 

50.490244a'3 - 120.995139r'3 + 12.624655a'V, (8) 



where 



and 



r' = i^-l. 



, RS 

a = 1. 

3V 



(9) 



(10) 



Values for R and S for the models studied in the work 
are given in Appendix A. 

The effect of varying c for a biaxial model is shown in 
Fig. [TJ where various models from Table |T] are plotted. 
As can be seen, the Vega EOS provides an excellent fit for 
all of the isotropic state points. This is also demonstrated 
in a plot of the self-dual models, where b = ^/ac (Fig. [2|) 
taken from Table [ill In figures [3] and [4] the isotropic 
equations of state are plotted for the 1 x 2.5 x 2.5, 1 x 1 x 
2.5 and for the 1x4x4, and 1x1x4 models. It can be 
seen that for these modest anisotropics the oblate/prolate 
curves are indeed almost coincident, as assumed in the 
simple EOS built using B2 only. The numerical results 
are presented in Tables IIIII and IIVI 



Bl = 18.3648 + 27.714434a' - 10.2046r'-h 11.142963a'' 



8.634491t'' 
17.190946a'2'' 



0.74674a'' 



28.279451a'r' 
r + 24.188979a'r'2 
9.455150t'3, 



(7) 



Isotropic-liquid crystal transition 

The theory behind the isotropic-nematic (I-N) transi- 
tion was first developed by Lars Onsager [43| and was 



TABLE I: Equations of state for biaxial (a = 1, a < b < c) hard ellipsoids. Error in y is 0(10 ''), error in Z is 0(10 ^) 
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FIG. 2: Plot of the simulation results for the self-dual 
molecules 1 x 1.25 x 1.5625, 1 x 2 x 4 and the isotropic state 
points for 1x3x9. The lines are the Vega equation of state. 
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FIG. 3: Plot of the results for the oblate 1 x 2.5 x 2.5 model 
(black squares) and the prolate 1 x 1 x 2.5 model (black circles) 
along with the Vega EOS (oblate: solid line, prolate: dashed 
line) 



TABLE II: Equations of state for biaxial self-dual (6 = ^/ac) 
hard ellipsoids. Error in y is O(10~^), error in Z is O(10~^) 



6= 1.25, c = 


1.5625 


b = 


2,c = 4 


6 = 3, c 


= 9 


* 

P 


y 


Z 


y 


Z 


y 


Z 


0.5 


0.141 


1.85 


0.126 


2.06 


0.104 


2.51 


1.0 


0.205 


2.55 


0.181 


2.87 


0.149 


3.50 


1.5 


0.247 


3.17 


0.218 


3.59 


0.181 


4.33 


2.0 


0.280 


3.73 


0.247 


4.23 


0.205 


5.10 


2.5 


0.303 


4.31 


0.271 


4.82 


0.228 


5.76 


3.0 


0.326 


4.81 


0.290 


5.40 


0.246 


6.38 


3.5 


0.344 


5.32 


0.305 


6.00 


0.264 


6.94 


4.0 


0.361 


5.79 


0.319 


6.55 


0.281 


7.44 


4.5 


0.375 


6.28 


0.335 


7.01 


0.301 


7.82 


5.0 


0.385 


6.79 


0.345 


7.57 


0.315 


8.29 


5.5 


0.389 


7.22 


0.359 


8.01 


0.332 


8.67 


6.0 


0.406 


7.72 


0.366 


8.56 


0.344 


9.12 


6.5 


0.420 


8.09 


0.376 


9.04 


0.356 


9.54 


7.0 


0.428 


8.55 


0.383 


9.56 


0.364 


10.06 


7.5 


0.437 


8.97 


0.394 


9.96 


0.376 


10.41 


8.0 


0.441 


9.49 


0.400 


10.45 


0.384 


10.88 


8.5 


0.451 


9.86 


0.408 


10.90 


0.395 


11.26 


9.0 


0.458 


10.28 


0.412 


11.42 


0.406 


11.58 


9.5 


0.463 


10.72 


0.419 


11.86 


0.410 


12.11 



studied for solutions of hard ellipsoids by Akira Isihara 
ill. Frenkel and Mulder found an I-N+ transition 

for a system of « 90— 108 prolate ellipsoids of c/a> 2.75. 
This result was called into question by Zarragoicoechea et 
al. [i^ , suggesting system size effects play an important 
role in locating the I-N transition. However, a later work 
by Allen and Mason [i^ confirmed the ncmatic phase for 
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FIG. 4: Plot of the results for the oblate 1x4x4 model (black 
squares) and the prolate 1x1x4 model (black circles) along 
with the Vega EOS (oblate: solid line, prolate: dashed line) 



c/a = 3.0 at densities of p/pcp ~ 0.73 — 0.75. In this work 
no such transition is found until c/a = 6. It is interest- 
ing to compare this to linear tangent hard spheres, where 
nematic phases were found for to = 6 and a smectic A 
phase for m = 5 where to is the number of monomer 
units [i^l (Note that smectic phases are not observed for 
hard ellipsoids). 

The appearance of the nematic phase for the 1x1x8 
(Fig. [7]) model was considerably 'delayed', and no ne- 
matic phase formed for compression runs (150k equilibra- 
tion followed by 150k production) of the prolate 1 x 1 x 10 
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TABLE III: Equations of state for oblate (a = 1, 6 = c) hard ellipsoids. Error in y is 0(10 error in Z is 0(10 ^) 
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TABLE IV: Equations of state for prolate (o = 6 = 1) hard ellipsoids. Error in y is 0(10 error in Z is 0(10 ^) 
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4.5 


0.353 


6.66 


0.327 


7.19 


0.314 


7.49 


0.304 


7.74 


0.285 


8.25 


0.269 


8.74 


5.0 


0.367 


7.12 


0.338 


7.73 


0.325 


8.03 


0.311 


8.40 


0.298 


8.77 


0.285 


9.15 


5.5 


0.377 


7.62 


0.348 


8.26 


0.332 


8.67 


0.326 


8.82 


0.311 


9.24 


0.297 


9.68 


6.0 


0.384 


8.17 


0.359 


8.74 


0.346 


9.06 


0.334 


9.39 


0.321 


9.76 


0.308 


10.19 


6.5 


0.396 


8.58 


0.367 


9.26 


0.352 


9.66 


0.349 


9.72 


0.330 


10.30 


0.317 


10.71 


7.0 


0.402 


9.10 


0.374 


9.79 


0.365 


10.03 


0.356 


10.27 


0.349 


10.49 


0.327 


11.17 


7.5 


0.411 


9.53 


0.385 


10.18 


0.373 


10.50 


0.377 


10.40 


0.356 


11.02 


0.334 


11.75 


8.0 


0.422 


9.92 


0.391 


10.70 


0.375 


11.16 


0.389 


10.76 


0.367 


11.40 


0.337 


12.40 


8.5 


0.428 


10.39 


0.395 


11.25 


0.388 


11.46 


0.402 


11.06 


0.376 


11.83 


0.345 


12.86 


9.0 


0.432 


10.90 


0.404 


11.65 


0.394 


11.94 


0.407 


11.55 


0.385 


12.21 


0.356 


13.20 


9.5 


0.440 


11.29 


0.413 


12.01 


0.403 


12.32 


0.417 


11.90 


0.389 


12.75 


0.367 


13.54 
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FIG. 5: Plot of the results for the oblate 1x5x5 model 
(isotropic: black squares, discotic: open squares) and the pro- 
late 1x1x5 model (black circles) along with the Vega EOS 
(oblate: solid line, prolate: dashed line) 



FIG. 7: Plot of isotropic the results for the oblate 1x8x8 
model (isotropic: black squares, discotic: open squares) and 
the prolate 1x1x8 model (isotropic: black circles, nematic: 
open circles) along with the Vega EOS (oblate: solid line, 
prolate: dashed line) 
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FIG. 6: Plot of isotropic the results for the oblate 1x6x6 
model (isotropic: black squares, discotic: open squares) and 
the prolate 1x1x6 model (isotropic: black circles, nematic: 
open circles) along with the Vega EOS (oblate: solid line, 
prolate: dashed line) 
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FIG. 8: Plot of the isotropic results for the oblate 1 x 10 x 10 
model (isotropic: black squares, discotic: open squares) and 
the prolate 1 x 1 x 10 model (black circles) along with the 
Vega EOS (oblate: solid line, prolate: dashed line) 



model (Fig. [8]). Given that it is fully expected to see 
phases for this model this indicates that so called 'jam- 
ming' is a real problem for especially elongated prolate 
systems. The system finds itself grid-locked, and is un- 
able to reorientate into the energetically more favourable 
nematic phase within the time scale of the simulation. 
The equation of state of the glassy state (Fig. [8]) can be 
seen to be intermediate between the higher compressibil- 
ity factor of the isotropic phase, indicated by the Vega 



EOS, and the much lower values for the discotic branch 
of the 1 X 10 X 10 model. On the other hand, discotic 
phases readily appeared for the 1x5x5 model (Fig. [5]) 
upwards. 

Samborski et al. 



48| placed the I-N transition at 



0.370 < y < 0.388 for the 1 x 1 x 5 model, 0.333 < 
y < 0.351 for the 1 x 5 x 5 model, 0.203 < y < 0.222 
for the 1 X 1 X 10 model and 0.185 < y < 0.203 for the 
1 X 10 X 10 model. It is interesting to note that I-N transi- 
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FIG. 9: Snap-shot of the 1 x 3 x 10 model in the biaxial 
phase (here at p* — 8.0). The ellipsoids are colour coded 
with respect to their orientations (colour online). 
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y 

FIG. 10: Plot of the equation of state for the biaxial 1 x 3 x 
10 model (isotropic: black squares, discotic: open squares, 
biaxial: black triangles) along with the Vega EOS for the 
isotropic branch. 



tion sets in at lower volume fractions for oblate ellipsoids 
than for prolate models that have the same value of c. 
The same is seen in this work for the '6' (Fig. [6]) and 
'8' (Fig. [7]) models. This is probably due to the fact 
that the compressibility factor of the oblate models, for 
a given volume fraction, is higher than that of its cor- 
responding prolate counterpart. This triggers the earlier 
onset of the I-N_ transition. 

In this study only one convincing biaxial phase was 
identified, that of the 1 x 3 x 10 model (sec Fig. ^ for 
a snapshot. Fig. [TO] for the EOS). The biaxial phase 



formed at p* = 8.0, however, at p* = 5.0 the system 
had formed a discotic phase. This indicates that the 
formation of a biaxial phase is a two-stage process; first 
orientating along the short a axis, later followed by the 
long c axis i.e. having the phase transitions isotropic - 
discotic - biaxial. This is in line with the observation 
that discotic phases form at lower volume fractions than 
nematic phases (the biaxial phase can be seen as being 
composed of both a discotic and a nematic phase). 

The positions of the isotropic liquid-liquid crystal tran- 
sitions found in this work are presented in Table [V] A 
study of the I-N transition for biaxial ellipsoids was un- 
dertaken by Tjipto-Margo and Evans [i^ confirming the 
finding of Gelbart and Barboy (H^l that biaxiality reduces 
the first order nature of the I-N transition (this result was 
also confirmed for hard sphero-platelets by Somoza and 
Tarazona [EH). It appears that the additional degree of 
freedom associated with biaxiality increases the disorder 
in the nematic phase with respect to a uniaxial (Dooh) 
model. At the 'self-dual' point (where a : b = b : c, i.e. 
b = \faE) the transition becomes second-order (s^ . The 
isotropic-biaxial nematic transition is said to occur for 
models that are close to this self-dual point [2^ . For the 
uniaxial 1 x IGx 10 model a considerable jump can be seen 
in the volume fractions associated with the formation of 
the discotic phase (Fig. [8]), one of the hall-marks of a 
first-order transition. Meanwhile, for the biaxial 1x3x10 
model no such jump is seen. 



CONCLUSIONS 

Various ellipsoidal models have been subjected to 
Monte Carlo compression runs. The Vega equation of 
state is seen to perform very well for the isotropic phases 
of all of the models considered. The elongated uniaxial 
models show indications of the first-order nature of the 
isotropic-ncmatic transition. However, the formation of 
orientationally ordered phases for very long prolate mod- 
els is severely hindered by the formation of a glass like 
state. The oblate models form discotic readily, at lower 
volume fractions than their prolate partners. The biaxial 
phase seems to form in a two-stage process, first form- 
ing a discotic phase, followed by the biaxial phase upon 
orientation of the long axes. 

The authors should like to thank C. Vega for the pro- 
vision of the Perram-Wertheim overlap algorithm and 
N. G. Almarza for useful discussions. This work was 
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0505/ESP/0299 - CSICQFT (MOSSNOHO) of the D. 
G. de Universidades e Investigacin del Comunidad de 
Madrid. One of the authors (C. M.) would like to thank 
the CSIC for the award of an I3P post-doctoral contract. 
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TABLE V: Location of the isotropic-nematic transition {S2 > 
0.4) for hard elhpsoids (Note: a = 1). Order parameters are 
given for p* = 9.5. 



TABLE VI: Values for R, S, V, a and B2/V for hard elhpsoids 
(Note: a = 1). 



b 


c 


P* 


y 


phase 


S2{c) 


52(a) 


prolate 




1 




2.5 






I 


< 0.1 


< 0.1 


1 




4 






I 


< 0.1 


< 0.1 


1 




5 






I 


0.13 


< 0.1 


1 




6 


7.0 


0.356 


N+ 


0.69 


0.20 


1 




8 


7.5 


0.356 


N+ 


0.55 


0.18 


1 




10 






I 


0.30 


0.12 


oblate 




2.5 




2.5 






I 


< 0.1 


< 0.1 


4 




4 






I 


0.12 


0.24 


5 




5 


7.0 


0.374 


N- 


0.23 


0.81 


6 




6 


5.5 


0.336 


N- 


0.25 


0.87 


8 




8 


3.0 


0.257 


iV_ 


0.26 


0.92 


10 




10 


2.5 


0.240 


N- 


0.27 


0.95 


biaxial prolate 


{b< , 


/ac) 




2 




5 






I 


0.10 


< 0.1 


2 




6 






I 


0.21 


0.15 


2 




8 


6.0 


0.342 


N+ 


0.85 


0.25 


3 




10 


8.0 


0.389 


B 


0.48 


0.76 


biaxial self-dual 




1.25 


1.5625 






I 


< 0.1 


< 0.1 


2 




4 






I 


< 0.1 


< 0.1 


3 




9 


4.5 


0.30 


N- 


0.34 


0.81 


biaxial oblate 


{b> , 


/ac) 




3 




6 


9.0 


0.402 


N- 


0.23 


0.46 


3 




8 


6.0 


0.338 


N- 


0.27 


0.78 


5 




8 


4.5 


0.313 


N- 


0.25 


0.90 


5 




10 


3.5 


0.269 


N- 


0.30 


0.92 


8 




10 


3.0 


0.249 


N- 


0.27 


0.95 



APPENDIX A 

In Table IVll we provide a table of R, S, V, a and B2/V 
for the various models studied in this work. The values 
for R and S are obtained by evaluating the expressions 
derived by Singh and Kumar (HsI . 54 1. Thus the mean 
radius of curvature is given by 



^-2 



1 + Cc Uc 



(11) 



and the surface area by 



S = 27ra^ 



1 + V^c{l + eb) \ -F{ip, fca) + E{ip, fca) 



where F{ip, k) is an elliptic integral of the first kind and 
E{ip^ k) is an elliptic integral of the second kind, with the 



k 

u 


c 


R 


5 


V 


a 


B2/V 


sphere 




-1 
t 


1 


1 


47r 


47r/3 


1 


4 


prolate 




1 


2.5 


1.5919 


26.1518 


10.472 


1.32516 


4.97549 


t 


4 


2.26639 


40.4975 


167r/3 


1.82597 


6.4779 


1 
t 


5 


2.73397 


50.1925 


207r/3 


2.184 


7.552 


1 
t 


6 


3.20942 


59.9386 


Stt 


2.55136 


8.65409 


1 
t 


8 


4.17441 


79.5147 


327r/3 


3.30174 


10.9052 


t 


10 


5.15042 


99.151 


407r/3 


4.06377 


13.1913 


oblate 




2.5 


2.5 


2.0811 


50.0111 


26.1799 


1.32516 


4.97549 


A 

4 


4 


3.22269 


113.921 


647r/3 


1.82597 


6.4779 


r 




5 


3.99419 


171.78 


1007r/3 


2.184 


7.552 


D 


6 


4.76976 


241.985 


487r 


2.55136 


8.65409 


Q 


8 


6.32758 


419.657 


2567r/3 


3.30174 


10.9052 


iU 


10 


7.89019 


647.22 


4007r/3 


4.06377 


13.1913 


biaxial 




1.25 


1.5625 


1.28328 


20.1576 


8.18123 


1.05395 


4.16185 


2 


4 


2.52566 


63.4766 


327r/3 


1.59473 


5.7842 


2 


5 


2.96925 


78.2743 


407r/3 


1.84951 


6.54853 


2 


6 


3.42527 


93.1895 


167r 


2.11675 


7.35026 


2 


8 


4.36059 


123.218 


647r/3 


2.67234 


9.01701 


3 


6 


3.70789 


129.13 


247r 


2.11675 


7.35026 


3 


8 


4.60996 


170.448 


327r 


2.60536 


8.81607 


3 


9 


5.07182 


191.203 


367r 


2.85815 


9.57446 


3 


10 


5.53883 


212.0 


407r 


3.11474 


10.3442 


5 


8 


5.22725 


268.73 


1607r/3 


2.7946 


9.38379 


5 


10 


6.10332 


333.946 


2007r/3 


3.24386 


10.7316 


8 


10 


7.1304 


521.211 


3207r/3 


3.69682 


12.0904 



amplitude being 



and the moduli 



and 



VJ = tan 



fci 



Ec - £6 



leb{l + ec) 



where the anisotropy parameters, et and Cc, are 



(12) and 



£6=1-1 - 1, 

. a , 



ec = 1- -1. 



(13) 
(14) 

(15) 

(16) 
(17) 
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The volume of the eUipsoid is given by the well known 

(18) 



V = ^a6c. 



Note the symmetry between the prolate and oblate ellip- 
soids for the values of a, and thus for B2- 



N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, 

A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 
(1953). 

B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 
1208 (1957). 

W. W. Wood, Los Alamos Scientific Laboratory Report 
LA-2827 (1963). 

W. W. Wood, J. Chem. Phys. 52, 729 (1970). 

J. Vieillard-Baron, J. Chem. Phys. 56, 4729 (1972). 

J. A. Cuesta and D. Frenkel, Phys. Rev. A 42, 2126 

(1990). 

D. Frenkel, B. M. Mulder, and J. P. McTague, Phys. Rev. 

Lett. 52, 287 (1984). 

B. Evans, Molec. Phys. 100, 199 (2002). 

D. Frenkel and B. M. Mulder, Molec. Phys. 55, 1171 

(1985). 

B. Mulder and D. Frenkel, Molec. Phys. 55, 1193 (1985). 
M. J. Preiser, Phys. Rev. Lett. 24, 1041 (1970). 
L. A. Madsen, T. J. Dingemans, M. Nakata, and E. T. 
Samulski, Phys. Rev. Lett. 92, 145505 (2004). 

B. R. Acharya, A. Primak, and S. Kumar, Phys. Rev. 
Lett. 92, 145506 (2004). 

J. W. Perram and E. Praestgaard, Molec. Phys. 63, 1103 
(1988). 

C. McBride and C. Vega, J. Chem. Phys. 117, 10370 
(2002). 

G.-W. Wu and R. J. Sadus, Fluid Phase Equilib. 194, 
227 (2002). 

M. Murat and Y. Kantor, Phys. Rev. E 74, 031124 
(2006). 

C. D. Michele, A. Scala, R. Schilling, and F. Sciortino, 
J. Chem. Phys. 124, 104509 (2006). 

A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Still- 
inger, R. Connelly, S. Torquato, and P. M. Chaikin, Sci- 
ence 303, 990 (2004). 

A. Donev, F. H. Stillinger, P. M. Chaikin, and 
S. Torquato, Phys. Rev. Lett. 92, 255506 (2004). 
A. Donev, S. Torquato, and F. H. Stillinger, J. Comput. 
Phys. 202, 737 (2005). 

A. Doneva, S. Torquatoa, and F. H. Stillinger, J. Com- 
put. Phys. 202, 765 (2005). 

M. G. Basavaraj, G. G. Fuller, J. Fransaer, and J. Ver- 
mant, Langmuir 22, 6605 (2006). 



[24] W. Man, A. Donev, F. H. Stillinger, M. T. Sullivan, 
W. B. Russel, D. Heeger, S. Inati, S. Torquato, and P. M. 
Chaikin, Phys. Rev. Lett. 94, 198001 (2005). 

[25] M. P. Allen, Liq. Cryst. 8, 499 (1990). 

[26] P. J. Camp and M. P. Allen, J. Chem. Phys. 106, 6681 
(1997). 

[27] W. Wang, J. Wang, and M.-S. Kim, Computer Aided 

Geometric Design 18, 531 (2001). 
[28] J. W. Perram and M. S. Wertheim, J. Comput. Phys. 58, 

409 (1985). 

[29] J. P. Straley, Phys. Rev. A 10, 1881 (1974). 
[30] R. Eppenga and D. Frenkel, Molec. Phys. 52, 1303 
(1984). 

[31] C. Zannoni, in The Molecular physics of liquid crystals, 

edited by G. R. Luckhurst and G. W. Gray (Academic 

Press, 1979), chap. 3, p. 51. 
[32] T. Boublik and I. Nezbeda, Coll. Czech. Chem. Commun. 

51, 2301 (1986). 
[33] A. Isihara, J. Chem. Phys. 18, 1446 (1950). 
[34] A. Isihara and T. Hayashida, J. Phys. Soc. Jpn. 6, 40 

(1951). 

[35] A. Isihara and T. Hayashida, J. Phys. Soc. Jpn. 6, 46 
(1951). 

[36] I. Nezbeda, Chem. Phys. Lett. 41, 55 (1976). 

[37] J. D. Parsons, Phys. Rev. A 19, 1225 (1979). 

[38] Y. Song and E. A. Mason, Phys. Rev. A 41, 3121 (1990). 

[39] M. Rigby, Molec. Phys. 66, 1261 (1989). 

[40] M. S. Wertheim, Molec. Phys. 99, 187 (2001). 

[41] M. Maeso and J. Solana, Molec. Phys. 79, 1365 (1993). 

[42] C. Vega, Molec. Phys. 92, 651 (1997). 

[43] L. Onsager, Ann. (N.Y.) Acad. Sci. 51, 627 (1949). 

[44] A. Isihara, J. Chem. Phys. 19, 1142 (1951). 

[45] G. Zarragoicochea, D. Levesque, and J. Weis, Molec. 

Phys. 75, 989 (1992). 
[46] M. P. Allen and C. P. Mason, Molec. Phys. 86, 467 

(1995) . 

[47] C. Vega, C. McBride, and L. G. MacDoweU, J. Chem. 

Phys. 115, 4203 (2001). 
[48] A. Samborski, G. T. Evans, C. P. Mason, and M. Allen, 

Molec. Phys. 81, 263 (1994). 
[49] B. Tjipto-Margo and G. T. Evans, J. Chem. Phys. 94, 

4546 (1991). 

[50] W. M. Gelbart and B. Barboy, Accounts of Chemical 

Research 13, 290 (1980). 
[51] A. M. Somoza and P. Tarazona, Molec. Phys. 75, 17 

(1992). 

[52] B. Mulder, Phys. Rev. A 39, 360 (1989). 

[53] G. S. Singh and B. Kumar, J. Chem. Phys. 105, 2429 

(1996) . 

[54] G. S. Singh and B. Kumar, Annals of Physics 294, 24 
(2001). 



